Patrick Metzger - Institut für Medizinische Bioinformatik und Systemmedizin
1 Einführung
Wir bitten euch hier, die grundlegenden Schritte einer RNA-Sequenzierungs-Analyse durchzuführen. Wir stellen euch bereits die auf das menschliche Referenzgenom ausgerichteten Reads zur Verfügung, d.h. euer Ausgangspunkt sind die “Raw Counts” für alle Gene jeder Probe. Zusammen mit den Counts stellen wir euch die Annotationstabelle zur Verfügung, die die Proben-IDs mit den tatsächlichen biologischen Bedingungen verknüpft. Der Datensatz stammt aus einem wissenschaftlichen Projekt, an dem wir hier in Freiburg arbeiten. Der Datensatz enthält Expressionsdaten eines Rektumadenokarzinoms, welches Varianten in BRAFD594G, KRASG12A und TP53R175H trägt. Mit dem RNA-Seq-Experiment soll untersucht werden, welchen Einfluss die Medikamente Sorafenib und Trametinib auf das Expressionsmuster der Tumorzellen haben. Um die Vergleiche vornehme zu können wurde ebenfalls eine Kontrollprobe mit Dimethylsulfoxid (DMSO) generiert. DMSO wurde gewählt, da die Medikamente darin gelöst sind, um sie den Zellen zuverabreichen. Wir betrachten im folgenden also die biologischen Bedingungen DMSO, Sorafenib (Sora) und Trametinib (Tram). Für alle Bedingungen liegen zwei biologische Replikate vor.
Das Ziel dieser Aufgabe ist es die signifikant veränderten Gene zwischen der Kontrolle (DMSO) und den beiden Behandlungen (Sora und Tram) zu finden. Des weiteren interssiert uns auch, ob es Unterschiede zwischen den Behandlungen gibt. Dabei hilft uns in beiden Fällen die Analyse der differentiell exprimierten Gene (DEG).
Im folgenden beschreiben wir die Aufgaben und geben euch weitere Hinweise und Erklärungen.
Wir werden in den nächsten Wochen schrittweise folgende Punkte bearbeiten:
Importieren der Count Daten und erstellen einer Count Matrix
Verknüpfung der Daten zu tatsächlichen biologischen Bedingungen
Erstellen einer Tabelle, die alle Gene und ihre zugehörigen Gen-Identifikatoren enthält.
Berechnung einer Hauptkomponentenanalyse (PCA), um einen ersten Überblick über den Datensatz zu erhalten
Die DEG-Analyse
Vergleich zwischen den Medikamenten
Gene-Set Enrichment Analyse (GSEA)
Visualisieren und Exportieren der Ergebnisse
Daraus resultieren folgende Aufgaben/ Fragen für euch. Wir haben diese auf die drei verbleibenden Wochen aufgeteilt. Weitere Details findet ihr fortlaufend in diesem Dokument.
Aufgaben/ Fragen Woche 4
Vorbereitung:
Importiert die Count-Daten und erstellt daraus eine Count-Matrix
Erstellt eine Tabelle, die alle gemessenen Gene und all ihre verschiedenen IDs enthält.
Ordnet den Proben ihre biologische Bedingung zu
Erstellt das DESeq2 Objekt und normalisiert die Daten
Erstellt zwei QC Abbildungen, einmal mit den Roh-Counts und einmal mit den normalisierten Count-Werten.
Aufgaben/ Fragen Woche 5
DEG-Analyse:
Berechnet eine Principal Component Analysis (PCA) und stellt diese grafisch dar.
Baut auf dem DESeq2 Objekt auf und führt die “Differentially Expressed Genes” (DEG) Analyse durch.
Bestimmt die DEGs für die beiden Vergleiche Sorafenib vs DMSO und Trametinib vs DMSO
Wie viele DEGs konntet ihr pro Vergleich identifizieren? (Cutoff padj < 0.05)
Wie verteilen sich die Anzahl der DEGs auf hoch- bzw. runter-reguliert? (log2FoldChange > bzw. < 0)
Stellt die Anzahlen grafisch z.B. als Barplot dar.
Baut die signifikanten DEGs als Übersichtstabellen ins HTML ein.
Stellt die Ergebnisse der DEG-Analyse mit Hilfe eines Volcano Plots dar.
Aufgabe/ Fragen Woche 6
Vergleich der beiden Medikamente:
Findet Gemeinsamkeiten und Unterschiede zwischen den beiden Medikamenten im Vergleich zu DMSO.
Wendet einen Exakten Test nach Fischer an um aus den Genen veränderte Gensets/ Signalwege abzuleiten. Diese Art der Analyse nennt man funktionelle Analyse oder auch “Gene-Set Enrichment Analysis (GSEA)”.
Stellt die Ergebnisse grafisch dar.
1.1 R Pakete
Verwendete R Paket:
Code
#if (!require("BiocManager", quietly = TRUE))# install.packages("BiocManager")#BiocManager::install(version = "3.21")library(DESeq2)library(openxlsx)library(pheatmap)library(org.Hs.eg.db)library(ggplot2)library(FactoMineR)library(ggrepel)library(apeglm)library(ggrepel)library(viridis)library(tidyverse)library(kableExtra)library(pheatmap)library(genefilter)library(EnhancedVolcano)library(UpSetR)library(msigdbr)library(doMC) # falls es hiermit Probleme gibt könnt ihr das Paket auch weglassenlibrary(foreach) # falls es hiermit Probleme gibt könnt ihr das Paket auch weglassenlibrary(BiocParallel) # falls es hiermit Probleme gibt könnt ihr das Paket auch weglassen# optional -> increase speed of analysis due to multicore processingregister(MulticoreParam(4)) # -> diese Zeile löschen, wenn es Probleme mit den Paketen doMC, foreach und/ oder BiocParallel gab
Hinweis: Falls diese nicht installiert sind, können sie mit dem Paketmanager von bioconductorBiocManager::install() oder install.packages() installiert werden.
Bei nicht UNIX basierten Betriebssystemen, z.B. Microsoft Windows, müsst ihr aufpassen, da dort der Backslash als Trenner verwendet wird. Ihr könnte aber, wenn ihr nicht sicher sein, auch die Funktion file.path(..., fsep = .Platform\$file.sep) vewenden, z.B.
Code
file.path("~","Documents","Lehre","BIDS","Bioinformatik und Systembiologie","2024","BIDS_RNA_Seq_2024","MIRACUM_BIDS_Bioinformatik_Systembiologie_RNA_Sequenzierung_Aufgabe_2024")
1.3 Funktionen
Wir stellen euch ein paar Funktionen zur Verfügung, die euch bei der Konvertierung der verschiedenen Gen-IDs unterstützen können. Alle Konvertierungen sind in der Hauptfunktion getGeneMat() zusammengefasst. Diese nimmt als Input die ENSEMBL Gen-IDs und konvertiert diese in die anderen und fasst alles zu einem “data.frame” zusammen.
2.1.1 Importieren der Rohdaten aus dem Alignment und der Quantifizierung
Das ausrichten der Sequenzierungsschnipsel (Alignment der Reads) wurde mit dem Progamm STAR gemacht. STAR bietet außerdem die Möglichkeit auch gleich die Qunatifizierung der Expression vorzunemhen. Dabei wurde der Parameter --quantMode GeneCount verwendet. Hierzu wurde eine sogenannte GTF/ GFF Annotationsdatei benutzt, welche die Information beinhaltet, welches Gene zu welchen chromosomalen Koordinaten gehört. Wir betrachten hier ein un-stranded-RNA-Sequenzierungs Experiment.
Die einzelnen .tab Dateien beinhalten die Counts pro Gene. Dabei gibt die 1. Spalte den ENSEMBL Gen-Identifier an, z.B. ENSG00000223972 und die 2. Spalte die entsprechenden un-stranded RNA-Seq Counts. Diese beiden Spalten brauchen wir im Folgenden. Die ersten vier Reihen geben ein paar Zusammenfassungsstatistiken über die Count Datei und werden nicht benötigt.
Erstellt die Count Matrix und die Genreferenztabelle mit allen Genen und den zugehörigen IDs. Stellt die Count-Matrix, die betrachteten Gene und die Annotation der Proben zu den biologischen Bedingungen innerhalb des HTML Dokuments dar.
Im nächsten Abschnitt wird das DESeq2-Modell erstellt. Damit wird die Grundlage für die spätere Analyse der differentiell exprimierten Gene (DEG) gelegt. Hierbei ist das Handbuch von DESeq2 sehr hilfreich. Wir importieren die Daten auf Grundlage der erstellten Count-Matrix, analog zu Count-Matrix-Import. Nach dem Import und dem Erstellen des DESeq2 Objektes müssen wir die Rohdaten einem Vorfilter Schritt unterziehen um sehr niedrige Counts zu entfernen. Hierbei is es hilfreich sich zu überlegen, wie viele biologische Replikate wir pro Bedingung haben. Idealerweise entfernt man alle Gene, die im vorliegenden Experiment in Summer über alle Bedingungen weniger als 5-10 Counts haben.
Wie viele Gene verlieren wir aufgrund von niedriger Expression?
sum(dim(count_matrix)[1]-sum(keep)) # discarding 39,575 records
[1] 39575
Code
dds <- dds[keep,]
2.1.3 Normalisierung und Differentielle Expression
Nachdem wir die schwach exprimierten Gene entfernt haben können wir die Normalisierung der Daten vornehmen. Dazu müssen wir die “size factors” und die “Dispersion” bestimmen, damit wir den tatsächlichen Signifikanztest (wir wollen den Wald- Test verwenden) anwenden können. Außerdem müssen wir die normalisierten Expressionswerte extrahieren, damit wir damit z.B. eine PCA erstellen können. Dazu verwenden wir die regularized logarithm (rlog) Transformation. Die Werte können dann mit der Funktion assay() extrahiert werden.
Zur besseren Interpretation und dem Verständnis der Normalisierung bietet es sich an die Count-Werte der einzelnen Gene als kombinierten Boxplot für die jeweiligen Bedingungen darzustellen; jeweils vor und nach der Normalisierung. Die Konvertierung der Count-Werte in log2(count + 1) hat sich in dieser Hinsicht bewährt.
Erstellt bitte einen Boxplot mit den log-transformierten Counts vor und einen Boxplot nach der Normalisierung.
Hinweis: Normalisierte Counts erhaltet ihr nach anwenden der DESeq() Funktion und der Extraktion der Werte mit counts(ddsObject, normalized = TRUE). Generell empfielt sich ggplot2 zum Zeichnen zu verweden, da es vielfälltige Möglichkeiten bietet und dadurch sehr ansehnliche Abbildungen erstellen kann. Um ein für ggplot2 passendes data.frame zu erstellen, können Pakete wie z.B. reshape2, dplyr, etc. hilfreich sein.
Code
## Roh-Counts extrahierenraw_counts <-counts(dds)# Log2 transformierenlog_raw <-log2(raw_counts +1)# Daten für ggplot umformendf_raw <-as.data.frame(log_raw) %>%rownames_to_column("gene") %>%pivot_longer(-gene, names_to ="sample", values_to ="log2count") %>%left_join(as.data.frame(colData(dds)) %>%rownames_to_column("sample"), by ="sample")# Boxplot erstellendf_raw %>%ggplot(aes(x = sample, y = log2count, fill = group)) +# why not group?geom_boxplot(outlier.size =0.3) +theme_bw() +labs(title ="Counts (trans) pro Bedingung",x ="Bedingung", y ="Counts (trans)") +theme(axis.text.x =element_text(angle =45, hjust =1))
Code
## Normalisierte Counts extrahierennorm_counts <-counts(dds, normalized =TRUE)# Log2 transformierenlog_norm <-log2(norm_counts +1)# Daten für ggplot umformendf_norm <-as.data.frame(log_norm) %>%rownames_to_column("gene") %>%pivot_longer(-gene, names_to ="sample", values_to ="log2normCount") %>%left_join(as.data.frame(colData(dds)) %>%rownames_to_column("sample"), by ="sample")# Boxplot erstellendf_norm %>%ggplot(aes(x = sample, y = log2normCount, fill = group)) +geom_boxplot(outlier.size =0.3) +theme_bw() +labs(title ="Counts (norm, trans) pro Bedingung",x ="Bedingung", y ="Counts (norm, trans)") +theme(axis.text.x =element_text(angle =45, hjust =1))
2.2 PCA - Principal Component Analysis
Bevor wir mit der tatsächlichen DEG-Analyse fortfahren berechnen wir eine PCA um einen ersten Eindruck von unseren Daten zu bekommen. Dazu verwenden wir die eben erstellten, normalisierten Expressionswerte (rlog-Werte). In der PCA sollen die Bedingungen als “individuals” betrachtet und die Gene als “variables”. Im Resultat soll jede Probe dargestellt und in seine Haupkomponenten zerlegt werden. Hierzu könnt ihr z.B. das Paket FactoMineR verwenden. Bitte erstellt einen sog. “Individuals Graph” mit den Proben und den ersten beiden Hauptkomponenten als Achsen. Wenn ihr das FactoMineR Paket verwendet findet ihr diese Infos unter pca$ind$coord. Um einen visuell ansprechenderen Abbildung zu erhalten würde ich empfehlen die Abbildung wieder mit dem Paket ggplot2 zu erstellen. Mit der Funktion ggsave() könnte ihr die mit ggplot2 erstellte Abbildung sehr einfach in ein geeignetes Format, z.B. PDF, png, etc. exportieren.
Erstellt eine PCA und stellt diese als Abbildung im HTML Dokument dar.
# Interpretation# ==============# PCA erklärt 87% der Gesamtvarianz in 2 Dimensionen# DMSO bildet einen Cluster, die Samples kann als (linear) homogen angenommen werden# Sorafenib bildet einen Cluster, die Samples können daher als (linear) homogen angenommen werden# Trametinib bildet zwei Cluster, die Samples sind daher nicht (linear) homogen### MDS: Multi-dimensional Scaling (linear relationship analysis)# Compute distance matrix from rlog-transformed datad <-dist(t(assay(rld))) # transpose because samples are columns# Perform classical MDS (multi-dimensional scaling)mds <-cmdscale(d, k =2) # k = number of dimensions (2D)# Create a data frame with sample informationmds_df <-as.data.frame(mds)mds_df$group <-colData(rld)$group # adjust 'group' to your grouping column# Plot with ggplot2ggplot(mds_df, aes(x = V1, y = V2, color = group)) +geom_point(size =3) +labs(x ="MDS1", y ="MDS2") +theme_minimal()
Code
# Interpretation# ==============# MDS zeigt ähnliches Bild wie PCA# Sorafenib bildet einen Cluster, die Samples können daher als (linear) homogen angenommen werden# DMSO und Trametinib bilden jeweils zwei Cluster, die Samples scheinen daher nicht-lineare Beziehungen zu haben, die von der MDS nicht erfasst werden### t-SNA: t-distributed stochastic neighbor embedding for non-linear relationship analysislibrary(Rtsne)# Extract assay data (genes x samples) and transpose to samples x genesmat <-t(assay(rld))# Optionally: reduce to top variable genes to speed up and denoisetop_var_genes <-head(order(apply(mat, 2, var), decreasing =TRUE), 500)mat_top <- mat[, top_var_genes]# Run t-SNE (you can adjust perplexity; try 5–30 depending on sample size)set.seed(42) # for reproducibilitytsne_out <-Rtsne(mat_top, dims =2, perplexity =1, verbose =TRUE)
Performing PCA
Read the 6 x 6 data matrix successfully!
Using no_dims = 2, perplexity = 1.000000, and theta = 0.500000
Computing input similarities...
Building tree...
Done in 0.00 seconds (sparsity = 0.611111)!
Learning embedding...
Iteration 50: error is 65.453220 (50 iterations in 0.00 seconds)
Iteration 100: error is 103.638598 (50 iterations in 0.00 seconds)
Iteration 150: error is 52.926096 (50 iterations in 0.00 seconds)
Iteration 200: error is 57.076918 (50 iterations in 0.00 seconds)
Iteration 250: error is 54.445751 (50 iterations in 0.00 seconds)
Iteration 300: error is 0.322697 (50 iterations in 0.00 seconds)
Iteration 350: error is 0.138587 (50 iterations in 0.00 seconds)
Iteration 400: error is 0.109798 (50 iterations in 0.00 seconds)
Iteration 450: error is 0.109646 (50 iterations in 0.00 seconds)
Iteration 500: error is 0.096316 (50 iterations in 0.00 seconds)
Iteration 550: error is 0.108245 (50 iterations in 0.00 seconds)
Iteration 600: error is 0.088462 (50 iterations in 0.00 seconds)
Iteration 650: error is 0.108734 (50 iterations in 0.00 seconds)
Iteration 700: error is 0.105920 (50 iterations in 0.00 seconds)
Iteration 750: error is 0.109760 (50 iterations in 0.00 seconds)
Iteration 800: error is 0.108264 (50 iterations in 0.00 seconds)
Iteration 850: error is 0.109467 (50 iterations in 0.00 seconds)
Iteration 900: error is 0.109448 (50 iterations in 0.00 seconds)
Iteration 950: error is 0.109438 (50 iterations in 0.00 seconds)
Iteration 1000: error is 0.109435 (50 iterations in 0.00 seconds)
Fitting performed in 0.00 seconds.
Code
# Combine with metadatatsne_df <-data.frame(TSNE1 = tsne_out$Y[,1],TSNE2 = tsne_out$Y[,2],group =colData(rld)$group # adjust if your metadata column is named differently)# Plot with ggplot2ggplot(tsne_df, aes(x = TSNE1, y = TSNE2, color = group)) +geom_point(size =3) +labs(title ="t-SNE of RNA-seq samples", x ="t-SNE 1", y ="t-SNE 2") +theme_minimal()
Code
# Interpretation# ==============# Ganz deutliche Cluster zu erkennen; DMSO und Tramatinib bilden jeweils einen Cluster, zusammen einen größeren. Daher kann davon ausgegangen, dass Trametinib keine über die Baseline hinausgehende nicht-lineare Beziehung aufweisen.# Sorafenib bildet einen Cluster deutlich abseits der Kontrollgruppe und von Trametinib, das kann als Indikator für einen tatsächlichen Unterschied gewertet werden
2.3 Differentielle Expressions Analyse
Um die DEG-Analyse durchzuführen müssen wir noch definieren an welchen tatsächlichen Vergleichen wir interessiert sind. Dazu können wir uns eine sog. Kontrast-Matrix erstellen. Folgende Vergleiche sind für uns interessant:
Sorafenib vs DMSO und
Trametinib vs DMSO,
da wir verstehen wollen, welchen Einfluss die beiden Inhibitoren auf unsere Zellen haben. Für beide Vergleiche bestimmen wir die DEGs und exportieren diese bis zu einem korrigierten (FDR) p-Wert < 0.05 in eine Tabelle. Hierzu verwenden wir die Funktion lfcShrink(). Damit die Tabellen für unsere Kollaborationspartner besser verständlich werden ist es essentiell, dass die Gene mit allen IDs, hautpsächlich aber dem Symbol, in den Ergebnissen enthalten sind.
Berechnet die DEGs für die beiden Vergleiche Sorafenib vs DMSO und Trametinib vs DMSO. Erstellt Tabellen, die die signifikanten Gene bis zu einem FDR-korrigierten p-Wert < 0.05 beinhalten. Achtet hier darauf auch das Symbol in der Ergebnistabelle zu haben. Ansonsten fällt die Zuordnung der Gene schwer. Zeigt diese im HTML Dokument. Exportiert die Ergebnisse der DEG-Analyse zusätzlich als Excel-Tabellen.
Bei der Darstellung der DEG-Ergebnisse innerhalb des HTML Dokuments müsst ihr nicht alle Spalten darstellen. Dies wird schnell sehr unübersichtlich. Ich würde sagen, dass ihr maximal drei Angaben braucht: Symbol, log2FoldChange und den korrigierten p-Wert. Bei der Darstellung innerhalb des HTMLs kann die Funktion kable() hilfreich sein.
Wie viele Gene sind pro Vergleich signifikant reguliert? Wie verteilt sich die Anzahl auf hoch- bzw. runter-regulierte Gene? Stellt diese Ergebnisse anschaulich dar.
Hinweis: Hierzu könnt ihr den log2FoldChange aus der DEG-Analyse verwenden.
Code
## ref:## https://www.youtube.com/watch?v=0b24mpzM_5M (theory)## https://www.youtube.com/watch?v=OzNzO8qwwp0 (R)if (!require("ashr", quietly =TRUE))install.packages("ashr")library(ashr)dds <-DESeq(dds)# Sorafenib vs DMSOres_sora <-lfcShrink(dds, contrast =c("group", "Sorafenib", "DMSO"), type ="ashr")# Trametinib vs DMSOres_tram <-lfcShrink(dds, contrast =c("group", "Trametinib", "DMSO"), type ="ashr")# Hilfsfunktion für Filter & Formatfilter_deg <-function(res) {as.data.frame(res) |>rownames_to_column("ENSEMBL") |>left_join(gene_reference[, c("ENSEMBL", "Symbol")], by ="ENSEMBL") |>filter(padj <0.05&!is.na(Symbol)) |>select(Symbol, log2FoldChange, padj) |>arrange(desc(abs(log2FoldChange))) # nach Effektgröße sortieren}deg_sora <-filter_deg(res_sora)deg_tram <-filter_deg(res_tram)## compare up and down-regulated genescount_regulation <-function(df) { df |>mutate(Regulation =case_when( log2FoldChange >0~"Upregulated", log2FoldChange <0~"Downregulated" )) |>count(Regulation)}sora_counts <-count_regulation(deg_sora)tram_counts <-count_regulation(deg_tram)# Visualizesora_counts$Comparison <-"Sorafenib vs DMSO"tram_counts$Comparison <-"Trametinib vs DMSO"all_counts <-bind_rows(sora_counts, tram_counts)ggplot(all_counts, aes(x = Comparison, y = n, fill = Regulation)) +geom_bar(stat ="identity", position ="dodge") +labs(title ="Differenziell regulierte Gene", y ="Anzahl Gene") +theme_minimal()
2.3.1 Identifizierte DEGs als Tabellen
Code
kable(head(deg_sora, 10), caption ="Top DEGs: Sorafenib vs DMSO")
Top DEGs: Sorafenib vs DMSO
Symbol
log2FoldChange
padj
TNS4
-4.276618
0e+00
FGF19
-3.548554
6e-07
SPRY4
-3.523185
0e+00
DUSP6
-3.409679
0e+00
PCK1
2.741600
0e+00
FOSL1
-2.659953
0e+00
PPP2R2C
-2.606917
0e+00
BNIP5
2.260004
0e+00
DUSP2
-2.197986
0e+00
KIAA0040
-2.185101
0e+00
Code
kable(head(deg_tram, 10), caption ="Top DEGs: Trametinib vs DMSO")
Zusätzlich zur tabellarischen Darstellung der DEG-Ergebnisse kann man diese auch als sog. “Volcano Plots” darstellen. Dabei werden alle Gene als Punkte mit ihrem Signifikanzwert (y-Achse) und dem log2FoldChange (x-Achse) dargestellt. Die signifikanten Gene werden dabei farblich hervorgehoben. Mit dem R Paket EnhancedVolcano haben wir gute Erfahrungen gemacht.
Erstellt jeweils einen Volcano Plot für die beiden Vergleiche.
Code
if (!requireNamespace("EnhancedVolcano", quietly =TRUE)) { BiocManager::install("EnhancedVolcano")}library(EnhancedVolcano)
Um die Unterschiede und Gemeinsamkeiten zwischen den beiden Medikamenten bzw. deren jeweiligen Vergleiche zu DMSO zu finden, können wir das R Paket UpSetR verwenden. Da wir beide Medikament zu DMSO verglichen haben, ist das Vorgehen sinnvoll und wir können direkt Unterschiede und Gemeinsamkeiten identifizieren. Wir nehmen den Vergleich auf Ebene der signifikant regulierten Gene, getrennt nach hoch- bzw. runter-Regulation, vor. Dazu legen wir uns zwei Listen an, analog zu Basic Usage. Eine für die hoch-regulierten und eine für die runter-regulierten Gene. Danach kann der Befehl upset(fromList(list)) verwendet werden. Um an die einzelnen Gene in den verschiedenen “Sets” zu kommen, stellen wir euch eine Funktion upSetSets(list) zur Verfügung. Als Input Parameter übergebt ihr dieser Funktion die Liste, die auch für upset() verwendet wird.
Führt den Vergleich der beiden Medikamente durch.
Code
if (!requireNamespace("UpSetR", quietly =TRUE)) {install.packages("UpSetR")}library(UpSetR)# function to obtain genes per setupSetSets <-function(sets){ list_names <-names(sets)attach(sets,warn.conflicts = F) res <-lapply(1:length(list_names),function(y){ combinations <-combn(list_names,y) res<-as.list(apply(combinations,2,function(x){if(length(x)==1){ p <-setdiff(get(x),unlist(sapply(setdiff(list_names,x),get))) }elseif(length(x) <length(list_names)){ p <-setdiff(Reduce(intersect,lapply(x,get)),Reduce(union,sapply(setdiff(list_names,x),get))) }else p <-Reduce(intersect,lapply(x,get))if(!identical(p,character(0))) pelseNA }))if(y==length(list_names)) { res[[1]] <-unlist(res); res<-res[1] }names(res) <-apply(combinations,2,paste,collapse="-") res }) result <-lapply(res, function(x) x[!is.na(x)]) result <-unlist(result, recursive = F) result <-lapply(result,function(x) data.frame(ID=x))return(result)detach(sets)}
2.4.1 Gemeinsamkeiten
2.4.1.1 Durch beide Medikamente hoch-regulierte Gene
Jetzt kennen wir die Gene, die entweder durch beiden Medikamente oder auch nur in dem einen oder dem anderen verändert sind. Aber was machen wir jetzt damit? Wir können uns z.B. der Ressource “Molecular Signatures Database” (MSigDB) bedienen.
“MSigDB is a resource of tens of thousands of annotated gene sets for use with GSEA (gene-set enrichment analysis) software”.
Die “gene-set enrichment” Analyse ist eine Berechnungsmethode, mit der festgestellt wird, ob ein a priori definierter Satz von Genen statistisch signifikante, übereinstimmende Unterschiede zwischen zwei biologischen Zuständen (z.B. Phänotypen) zeigt. Anders ausgedrückt können wir damit bestimmen, ob die Gene innerhlab der oben bestimmten “Sets” einen signifikante “Funktion” haben und daraus schlußfolgern, dass diese “Funktion” in unserem Experiment verändert ist. Für die Analyse verwenden wir den Exakten Test nach Fisher oder auch hypergeometrischer Test. Die Funktion hyperG() für den Test stellen wir euch zur Verfügung.
Code
hyperG <-function(geneSets,DEgenes,universe, cutoff=0.1, mincount=2, parallel=T, adj.P.Val = F,set.size =NULL){#' hyperG#' #' @description Calculates Fisher's Exact test with the specified genes and the supplied gene-sets.#'#' @param geneSets list. Gene-Set the calculation is based on, e.g. go.bp#' @param DEgenes character vector. Gene IDs used for testing. Same identifiers as used for the gene-sets, e.g. ENTREZ IDs.#' @param universe character vector. Universe gene IDs.#' @param cutoff numeric. Cutoff used to identify sig. pathways. Default: 0.1.#' @param mincount numeric. Consider only pathways which contain at least mincount genes. Default: 2#' @param parallel boolean. Use parallel calculation. Default: TRUE#' @param adj.P.Val boolean. Use adjusted p-value for significance filtering. Is always calculated.#' @param set.size vector. Min and max size of allowed gene-sets. Default min:10 genes and max:500 genes.#' #' @return the significant regualted pathways.#' @export#' @importFrom foreach, doMCrequire(foreach)require(doMC)if(parallel){registerDoMC(cores=detectCores()) cores=detectCores() }else{ cores=1 }if(!is.null(set.size)){print('Set Size Limits') idx <-lapply(geneSets,function(x){length(x) <= set.size[2] &length(x) >= set.size[1]}) geneSets <- geneSets[unlist(idx)] } l <-length(setdiff(universe,DEgenes)) DElen <-length(DEgenes) results <-mclapply(1:length(geneSets), function(i){ results <-matrix(data=NA,ncol=7,nrow =1)colnames(results) <-c('Term','Count','Size','p-value','adj.P.Val','odds ratio','GeneIDs') geneSet <-intersect(universe, geneSets[[i]]) e <-intersect(DEgenes,geneSet) a <-length(e) b <- DElen - a c <-length(geneSet) - a d <- l - c contigency.matrix <-cbind(c(a,b),c(c,d)) res <-fisher.test(contigency.matrix,alternative ='greater') results[1,'Term'] <-names(geneSets)[i] results[1,'Count'] <- a results[1,'Size'] <-length(geneSets[[i]]) results[1,'p-value'] <- res$p.value results[1,'odds ratio'] <- res$estimate[[1]]# find genes annotated in the consensus termif(a >0){ genes <-intersect(DEgenes,geneSet) eid <- genes eid <- eid[order(eid)] results[1,'GeneIDs'] <-paste(eid,collapse="|") }return(results) }, mc.cores=cores) results <-as.data.frame(do.call(rbind, results))for(i inc(2, 3, 4, 5)){ results[, i] <-as.numeric(as.character(results[, i])) }if(nrow(results) !=1){ results <- results[order(results[,'p-value'],decreasing =FALSE),] results[,'adj.P.Val'] <-p.adjust(results[,'p-value'], 'BH')if(adj.P.Val){ results <-as.data.frame(subset(results,results[,'adj.P.Val']<=cutoff)) }else{ results <-as.data.frame(subset(results,results[,'p-value']<=cutoff)) } results <-as.data.frame(subset(results,results[,'Count']>=mincount)) }else results <-as.data.frame(results)return(results)}
Außer der Funktion benötigen wir noch die entsprechenden Gensets/ Signalwegen, die uns interessieren. Diese können wir direkt über ein R Paket von MSigDB (msigdbr) beziehen. Auch hierzu stellen wir euch eine Funktion zur Verfügung, die das abrufen und erstellen der passenden Gensets vereinfacht. Wir laden uns damit die “Hallmark Gene Sets” direkt mit den Gen-Symbolen der Gene als IDs. Das geladene Genset ist eine Liste, die direkt als Input für die Funktion hyperG() verwendet werden können.
Um die Funktion hyperG() ausführen zu können brauchen wir jetzt noch ein sog. “universe”. Dieses beinhaltet alle Gene, die in unserem Experiment enthalten sind. Das “universe” ist ein Vektor, der alle Gene als Symbole beinhaltet. Es sollten keine NAs und/ oder Duplikate enthalten sein.
Erstellt das “universe” und führt die funktionelle Analye/ Gene-set Enrichment Analyse mit der hyperG Funktion und den Hallmark Signalwegen durch
Code
# Extrahiere die ENSEMBL-IDs aus dem dds-Objektensembl_ids <-rownames(dds)# Verknüpfe mit gene_reference, um die Symbolnamen zu erhaltenuniverse_df <-tibble(ENSEMBL = ensembl_ids) |>left_join(gene_reference[, c("ENSEMBL", "Symbol")], by ="ENSEMBL")# Erstelle das Universeuniverse <- universe_df$Symbol |>unique() |>na.omit()
2.4.1.3 Gemeinsamkeiten: Durch beide Medikamente hoch-regulierte Hallmark Gensets
## Interpretation (KI-unterstützt)## ==============## Die Behandlung mit Sorafenib führt zu einer signifikanten Hochregulation von Genen, die mit dem Zelltod (Apoptosis, p53), zellulärem Stress (zB Uv-Antwort), Immunantwort (Interferon Gamma) und Stoffwechselprozessen (Adipogenese, Häm- und Peroxisomenstoffwechsel) assoziiert sind. Dies unterstützt die bekannte cytotoxische und wachstumshemmende Wirkung von Sorafenib auf zellulärer Ebene.enrichment_up_tram <-hyperG(hallmark, up_tram, universe, cutoff =0.05, adj.P.Val =TRUE) # empty???top_enrichment_plot(enrichment_up_tram, "Top Enriched Hallmark Pathways (Trametinib UP)")
Code
## Interpretation## ==============## ?
2.4.1.4 Gemeinsamkeiten: Durch beide Medikamente runter-regulierte Hallmark Gensets
## Interpretation (KI-unterstützt)## ==============## Sorafenib führt zur signifikanten Herunterregulierung von Genen, die an Zellwachstum (MYC, mTORC1, G2M), Zellzyklusprogression, Entzündungsreaktionen (TNFα/NFκB) und Stressantworten beteiligt sind. Diese Ergebnisse sprechen für eine starke antiproliferative und antiinflammatorische Wirkung des Medikaments, ergänzt durch eine Hemmung zellulärer Überlebenssignale.enrichment_down_tram <-hyperG(hallmark, down_tram, universe, cutoff =0.05, adj.P.Val =TRUE)top_enrichment_plot(enrichment_down_tram, "Top Enriched Hallmark Pathways (Trametinib DOWN)")
Code
## Interpretation (KI-unterstützt)## ==============## Trametinib führt zu einer ausgeprägten Unterdrückung entzündungsbezogener Signalwege (TNFα/NFκB, IL-2, IL-6, STAT3), Zellüberlebenssignale sowie Komponenten der MAPK-Kaskade (KRAS). Diese Ergebnisse sind kompatibel mit dem Wirkmechanismus von Trametinib als MEK-Inhibitor und deuten auf eine kombinierte hemmende Wirkung auf Proliferation, Entzündung und Stressantwort hin.
